Reducing orbital eccentricity of precessing black-hole binaries 
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Building initial conditions for generic binary black-hole evolutions which are not affected by initial 
spurious eccentricity remains a challenge for numerical-relativity simulations. This problem can be 
overcome by applying an eccentricity-removal procedure which consists of evolving the binary black 
hole for a couple of orbits, estimating the resulting eccentricity, and then restarting the simulation 
with corrected initial conditions. The presence of spins can complicate this procedure. As predicted 
by post-Newtonian theory, spin-spin interactions and precession prevent the binary from moving 
along an adiabatic sequence of spherical orbits, inducing oscillations in the radial separation and in 
the orbital frequency. For single-spin binary black holes these oscillations are a direct consequence 
of monopole-quadrupole interactions. However, spin-induced oscillations occur at approximately 
twice the orbital frequency, and therefore can be distinguished and disentangled from the initial 
spurious eccentricity which occurs at approximately the orbital frequency. Taking this into account, 
we develop a new eccentricity-removal procedure based on the derivative of the orbital frequency 
and find that it is rather successful in reducing the eccentricity measured in the orbital frequency 
to values less than 10~* when moderate spins are present. We test this new procedure using 
numerical-relativity simulations of binary black holes with mass ratios 1.5 and 3, spin magnitude 
0.5, and various spin orientations. The numerical simulations exhibit spin-induced oscillations in 
the dynamics at approximately twice the orbital frequency. Oscillations of similar frequency are 
also visible in the gravitational-wave phase and frequency of the dominant I = 2, m = 2 mode. 

PACS numbers: 04.25.D-, 04.25. dg, 04.25.Nx, 04.30.-w 
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I. INTRODUCTION 

Over the last few years, numerical simulations of bi- 
nary black holes have improved tremendously (see e.g., 
the recent reviews [H-Q ) . These simulations are now used 
to aid data analysts for gravitational-wave detectors in 
the construction of analytical templates , and in test- 
ing the efficienc y of detector pipelines by injecting numer- 
ical waveforms [lG |. 

During the gravitational-radiation driven inspiral of a 
binary black hole, the orbital eccentricity decreases very 
quickly (ill , [l^ . For binary black holes formed from 
binarystellar evolution [l3| (instead of dynamical cap- 
ture [l3, EBDj the orbital eccentricity is expected to be 
essentially zero by the time the binary enters the fre- 
quency band of ground-based gravitational-wave detec- 
tors. Therefore, it is important that numerical simula- 
tions can be done for very low eccentricity binaries. 

Performing black-hole simulations with very small or- 
bital eccentricity is not easy for several reasons. Orbital 
parameters that result in vanishing eccentricity arc only 
known approximately through post-Newtonian (PN) the- 
ory . The translation of orbital parameters from PN 
theory into a complete binary black-hole initial-data set 
is ambiguous, because of differing coordinate systems and 
effects arising from solving the non-linear Einstein con- 
And finally, early in a numerical 



17 



straint equations 
evolution each black hole relaxes toward a steady state, 
affecting the black- hole masses, spins [l"8l - [20j . and orbital 
parameters. 



The complete evolution of a binary black hole is de- 
termined by its initial data. Therefore, control of orbital 
eccentricity has to be addressed in the construction of 
the initial data. The first approaches to construct low- 
eccentricity initial data were based on the assumption of 
circular orbits with the orbital fre quen cy determined by 
the effective potential method [2l| - |23| and the "Komar 
mass" ansatz [24[-p 8l. Both methods were shown to give 
similar results (28l| . Reference [29j presented techniques 
to measure eccentricity based on initial data alone. When 
binary black- hole evolutions became possible (sol - fs^ . it 
was realized that initial data constructed using the as- 
sumption of circular orbits resulted in a spurious orbital 
eccentricity of order one percent [33l - l35| . primarily due 
to neglecting the radial inspiral velocity, and due to the 
initial relaxation of the black holes. 

Two techniques are in use to achieve an orbital ec- 
centricity smaller than what can be achieved with quasi- 
circular initial data. One approach [s^ evolves PN equa- 
tions for the trajectories of the centers of the black holes. 
This subsidiary evolution of ordinary differential equa- 
tions is started at large initial separation, so that any 
spurious eccentricity due to the initial conditions dies 
out and the binary settles down to an inspiraling orbit 
with non-zero radial velocity. At the desired separation, 
the subsidiary evolution is stopped, the positions and 
velocities of the particles are read off, and are used as 
initial data parameters for the construction of the initial 
data set for the subsequent numerical evolution. This 
approach reduces eccentricity to about 0.002 for equal- 
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mass, non-spinning binaries, but is less successful for un- 
equal masses or high spins |37| . 

The second approach, proposed in Ref. (s^l and re- 
fined in Ref. |38| . performs an iterative procedure (see 
also Ref. [S^ ) . One begins with initial data with reason- 
ably low eccentricity, e.g., quasi-equilibrium initial data 
or initial data utilizing PN information. One evolves this 
initial data for about two to three orbits, analyzes the 
orbit, and computes an improved initial data set with 
(hopefully) lower eccentricity. This procedure can be re- 
peated until the desired degree of eccentricity is obtained. 

In past applications, eccentricity removal was based on 
the behavior of the proper separation s between the black 
hole apparent horizons. For non-spinning black hole bi- 
naries [40| or binaries with spins parallel to the orbital 
angular momentum pol . l4l| . this works very well, and 
the eccentricity drops by about an order of magnitude 
with each iteration. However, when one applies this ec- 
centricity removal procedure to precessing binaries, one 
encounters the difficulties illustrated in Fig. [T] At high 
eccentricity e > 0.01, s (where we indicate with a dot 
the derivative with respect to time) shows the expected 
oscillations with a period somewhat longer than the or- 
bital pcriod0 At sufficiently small eccentricity, however, 
the proper separation s and the radial velocity s exhibit 
oscillations at twice the orbital frequency (or one-half the 
orbital period). This frequency is distinct from the fre- 
quency of oscillations caused by orbital eccentricity, and 
its presence makes it very hard to further reduce eccen- 
tricity based on an analysis of the proper separation s. 

In this paper, we investigate these oscillations at twice 
the orbital frequency, and develop techniques for eccen- 
tricity removal for precessing binaries that mitigate the 
issues illustrated in Fig. [TJ We can understand these os- 
cillations based on PN calculations. In fact, as it was 
shown in Refs. [43l - l45| . spin-spin interactions in the dy- 
namics and spin precession, can introduce oscillations in 
the orbital separation and orbital frequency that pre- 
vent the binary from moving along a sequence of spher- 
ical orbits. Moreover, for single-spin binary black holes, 
the presence of spin-induced oscillations in the dynamics 
is a direct consequence of monopole-quadrupole interac- 
tions, that is of interactions of the form mi S'|/m2 and 
TO2 Sf /nil ■ It turns out that the amplitude of the spin- 
induced oscillations in the orbital frequency is half the 
amplitude of the oscillations in the coordinate separation. 
Therefore, we propose to base the iterative eccentricity 
removal on the orbital frequency and its time-derivative. 
We develop the relevant updating formulae for iterative 
eccentricity removal based on the (time-derivative) of the 
orbital frequency, and demonstrate with fully numerical 
simulations that iterative eccentricity removal can pro- 
ceed to much smaller eccentricities e < 10^* that are 
measured in the orbital frequency. 
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FIG. 1. Eccentricity removal based on proper separa- 
tion applied to a precessing binary black hole: The horizon- 
tal axis represents time in units of the initial orbital period 
P = 442M. For eccentricity e > 0.01, oscillations due to 
orbital eccentricity with period ~P dominate, and eccentric- 
ity removal is effective. For e ~ 0.01 oscillations at one-half 
the orbital period become apparent, spoiling further eccen- 
tricity removal based on s. In this example, the mass-ratio is 
mi/m2 = 1.5, the larger black hole carries spin — 0.5mf 
with initial direction tangent to the orbital plane, and the 
smaller black hole has zero spin. 



^ The period of radial oscillations exceeds the orbital period be- 
cause of periastron advance [43l . 



We also find that PN theory predicts spin-induced os- 
cillations in the separation with much smaller amplitude 
than those visible in Fig. [T] Figure [1] utilizes the proper 
separation s between the apparent horizons along a line 
element joining their centers. We find that use of the co- 
ordinate separation between the centers of the apparent 
horizons instead results in much smaller oscillations. Fi- 
nally, we find that the spin-induced oscillations arc also 
present in the gravitational-wave frequency and phase, 
and are qualitatively reproduced by the simple PN model 
used here. 

This paper is organized as follows. In Sec. [Ill we work 
out the spin-induced oscillations in the radial separation 
and orbital frequency for a PN model utilizing the Taylor- 
expanded PN Hamiltonian with only the lowest order PN 
terms responsible for the spin-induced oscillations. We 
also compare the obtained analytic formulae with nu- 
merical solutions of the ordinary differential equations 
describing a PN binary. In Sec. IIIIl we present the new 
method to reduce orbital eccentricity in presence of spins. 
In Sec. llVi we apply our improved eccentricity-removal 
procedure to fully general-relativistic simulations of sin- 
gle and double spin binary black holes, and compare 
the results with the earlier eccentricity-removal proce- 
dure based on the proper horizon separation. We also in- 
vestigate the presence of spin-induced oscillations in the 
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gauge-invariant gravitational-wave phase and frequency 
of a single-spin binary black hole. Finally, in Sec. |V] we 
summarize our main conclusions. 



II. ECCENTRICITY AND SPIN-INDUCED 
OSCILLATIONS 

As mentioned in the introduction, spin-spin effects and 
precession can induce oscillations in the orbital radial 
separation and frequency preventing the binary black 
holes from moving along an adiabatic sequence of spher- 
ical orbits. This result can be obtained in a straightfor- 
ward manner in PN theory (43l - l45l |. Here we re- analyse 
it in some detail using the PN Hamiltonian formalism. 

As we shall see, when spins are not aligned with the 
initial direction of the orbital angular momentum, the 
spin-induced oscillations are unavoidable and their im- 
portance increases at smaller distances, since at leading 
order spin-spin interactions scale as where r is the 

binary separation. Thus, if we were to start the binary 
black-hole evolution at large separation with some ini- 
tial orbital eccentricity, we expect that, by the time the 
binary reaches smaller separations, only spin-induced os- 
cillations would be left. However, numerical-relativity 
simulations start the evolution at a separation where 
the initial orbital eccentricity is not negligible. As we 
shall discuss, there exists an efficient way to distinguish 
and disentangle the initial orbital eccentricity from the 
spin-induced oscillations, namely the typical frequency 
at which these two effects occur. 

Henceforth we use natural units G = c = 1. 



A. Eccentricity in Newtonian dynamics 

Here wc briefly review some useful formulae of Newto- 
nian dynamics of eccentric orbits that we shall use below. 

In the ccnter-of-mass frame, the two-body problem 
reduces to a one-body problem for a particle of re- 
duced mass jjL = nil 1712/ M, subject to the acceleration 
r = —M/r^ r, where M = mi + TO2 is the total mass of 
the binary. In the Keplerian parametrization, a Newto- 
nian orbit of eccentricity e can be described in terms of 
the eccentric anomaly u (see, e.g., Ref. [46j ) 

u{t) — esmu{t) ^ Clt, (1) 
where Cl = \J~Mja^ , a being the semi-major axis, so that 



r{t) 



[1 — e cos ■ti(i)] 



(2) 



where r = a(l — e^). In the limit of small e we can 
approximate r[t) using only the first harmonic, that is 



"[t) =f [l-e cos(m)] +0{e^). 



(3) 



In fact, although the frequency spectrum of r{t) con- 
tains all harmonics in fJ, a Fourier analysis of r{t) shows 



that harmonics beyond the first one are quite suppressed 
in presence of a small eccentricity (see, e.g., Ref. H^). 
In the Keplerian parametrization the orbital frequency 
reads 



n{t) 



(4) 



[1 — e cos u{t)\ 

and in the limit of small e we find 

Vt{t)=Cl[l + 2ecos{nt)\+0{e^). (5) 

So, in Newtonian dynamics, whenever we have an eccen- 
tricity in rit), we expect oscillations of amplitude 2e f2 at 
the frequency Cl in ^l{t). 



B. Two-body dynamics for spinning black holes in 
PN theory 

We consider a binary composed of two black holes with 
masses mi and m2 and spins Si and S2. The binary 
dynamics can be described using the spinning Taylor- 
expanded PN Hamiltonian. In the center-of-mass frame, 
the Hamiltonian depends on the canonical variables (r, p) 
which describe the motion of a particle of reduced mass 
/i. and on the spins Si and S2. 

For the purposes of our analysis, it is sufficient to 
restrict the discussion to the Newtonian Hamiltonian, 
Hjsicwt, and include only the leading 1.5PN spin-orbit 
(SO) interaction [47j. Hsp, and the leading 2PN spin- 
spin (SS) interaction [4^, Hgs , where the SS interaction 
includes spin- induced monopole-quadrupole terms [49j . 
The Hamiltonian reads 



H — H{r, p; Si, S2) — i?Ncwt + ^^so + Hi 



SS, 



where 



-ffNowt — 



p" /.iM 



2fi r ' 



^^so = ^Soff • L, 



Hi 



SS 



with L = r X p, n = r/|r| and 
3to2 



3 {So -hf -Si 



Soff ^ 

So 



1 



4m 1 



Si+ 1 



3mi 

4m2 



1 + — I Si+ 1 + — ) S2. 

mi J \ m2, 



(6) 

(7) 

(8) 
(9) 

(10) 
(11) 



For reference, we point out that Sq can be rewritten in 
dimensionless form as 



So _ mi Si m2 S2 



(12) 
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The Hamilton equations of motion are given by 

OPi 



Pi 



(13) 
(14) 



where Fi is the radiation-reaction force. Here we foUow 
Ref. [13] and express Fi in terms of the Newtonian en- 
ergy flux J^n = 32/zV(5M2)?;io [seeEqs. (3.15), (3.27) in 
Ref. [H^l)] where for quasi-circular orbits v = (il/il)^/'^. 
Equations (fTS]) . (fT4|l must be supplemented with the spin 
precession equations 



dFl 



dS{ 
dH 



(15) 
(16) 



where e^^^ is the Levi-Civita symbol in flat spacetime. 
The non-spinning conservative part of the dynamics to- 
gether with the lowest-order SO interactions allows the 
existence of spherical orbits [r(t)=const.] [i^. In fact, 
if we consider H = i?Ncwt + ^^soi it is straightforward 
to show that the Hamiltonian is a spherically symmet- 
ric function that depends only on the radial separation 
and its conjugated momentum, i.e., H = H(r,pr). This 
happens because the other degrees of freedom are con- 
strained by the constants of motion: and Seff ■ L. 
More explicitly 



Imposing that at t = 0, r ~ vq = const., we have 



'dH' 




'Pr_' 


dpr 








= 



(17) 



(18) 



and to have a stable spherical orbit we have to require 
also that [pr]o = 0, hence 



[Pr]a 



dH 

dr 

ill 



{r,Pr = 0) 



Soft • L 



= 0. 



(19) 



Choosing properly and Soff • L to satisfy Eq. (fT9)) . we 
obtain spherical orbits. However, once SS interactions 
are included, this is no longer true, since and SeS-L are 
no longer constants of motion. Therefore, we must expect 
oscillations induced by SS terms in the radial separation 
and orbital frequency about their average values. 



C. Oscillations induced by leading SS interactions: 
conservative dynamics 

In this section, we investigate the behavior of the ra- 
dial separation r and of the orbital frequency fi at 2PN 



order for the conservative non-spinning dynamics. While 
doing so, we will also assume a negligible precession of 
the spins and of the orbital plane, since it takes place on 
a longer timescale than the effects we are interested in. 
Henceforth, we follow the method outlined in Appendix 
B of Ref. [111. 

As a first step, we restrict ourselves to the case in which 
radiation-reaction is not present (i.e., Fi = 0). As dis- 
cussed earlier, the presence of SS terras prevents r and 
from being constant. Thus, we write [45| 



(20) 



where the bar stands for time-average (...) over one or- 
bital period; hence, by definition, {Sr{t)) ~ {5^{t)) = 0. 
Our goal is to determine the equations that the oscil- 
lations 5r{t) and Sn{t) must obey at 2PN order. For 
convenience, we decompose vectors with respect to the 
triad defined by 



r X r 



(21) 



This triad is such that fi and A are in the instantaneous 
orbital plane, while Ln is orthogonal to it. In the instan- 
taneous orbital plane, we have the velocity 

V = r = rh + QrX, (22) 

and the acceleration 

a = arad n -|- Otan A -(- a_LLN , (23) 

with 

Orad = n • f = f- rl7^ (24) 



otan = A • i= = --^ (r^ n) , 



and 



a I = Lm ■ r = rfl I \ ■ 



(25) 



(26) 



For future reference, we define the projection of So on 
the instantaneous orbital plane 



So — So • Lj, 



(27) 



Note that Eq. (p2|) implicitly defines fl. We have = 
(r • A)/r. Wc need the acceleration r so we take a time 
derivative of Eq. and substitute Eq. into that. 
We note that the Newtonian orbital angular momentum 
can be written as 

Ljsi = tinr^L^, (28) 

while the orbital angular momentum L = r x p can be 
obtained from the Hamilton equation (jl3|) . that is to say 



P 



(Seff X r) 



(29) 
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so that 



2m 

r 



n(S 



off ■ n) 



(30) 



Since wc want to work consistently at 2PN order, we 
replace L with Eq. pop whenever it shows up in our cal- 
culations and drop higher PN terms. It is then straight- 
forward to compute the radial component of the acceler- 
ation 



Orad 



^ < 1 



2Mh 



3 (So 



■(Soff • Ln) 



(31) 



Since the leading-order spin acceleration is of 1.5PN 
order, we assume that the radial oscillations scale at least 
as r = 0(3)0 Hence, when computing the tangential 
component of the acceleration, at 2PN accuracy, we set 
r = (v = r f2 A) and pr = m every term coming 
from SO or SS interactions. Moreover, we also neglect 
any term depending on the time derivative of the spin in 
atan, since it is of higher PN order. This means that we 
are implicitly assuming that the spins are constant. At 
2PN order, we are left with 

ata„ = --r^(So-ft)(So-A). (32) 

Combining Eq. (^5]) with Eq. ([5^ . we solve for O by 
integrating Eq. ([5^ . In doing that, we assume that r 
and O in the right-hand side of Eq. ([32]) are constants, 
as their time derivatives are at least 0(3), and also the 
spins are constants, that is they do not precess. Thus, 
under those assumptions, the time evolution of the triad 
{n. A, Ln} is such that n and A swipe the orbital plane 
at a frequency Q while Ln stays fixed: 

h{t) = cos{nt) no + sm{nt) Aq, (33) 

X{t) = - sm{nt) ho + cos(m) Ao, (34) 

Ln(0 = Lno, (35) 

where iio = n(0), Ao = A(0) and Lno = Ln(0), and also 



n = r^A 



A = r^n. 



(36) 



Moreover, since we assume that the spins remain con- 
stant, we formally set Si(t) = Si(0) and 82(1) = 82(0), 
so in what follows So = So(0) and Soff = Sefi(O). By 
inserting Eqs. ((20)l into Eqs. (|3T|) and ([32)) . one obtains 
a pair of coupled differential equations: 



6r 



2n 

— (Soff • Ln) 



3 (So 



o2 
^0 



M 



(37) 



We denote the nPN order as 0{2n) 



and 



2nfSr 



f2 sn 



2M n 1 



r(So-fi)' 



(38) 



Here, k is an integration constant and, again, in the right- 
hand-side of the above equations we keep only terms 
through 2PN order. To fix k we time-average the above 
equations. We have 



((So • m) 



obtaining 



1 



NO ^NO 







Sq — (So • Lno) 



AM Vt f 3 



^0 ^ (^0 • Lno)" 



(39) 



(40) 



Taking the time average of Eq. ([37)) . we derive the 
following modified version of Kepler's law relating f and 



-,_M 2n . 

" — ^ — rr^off • -Lino + 



4Aff5 



3(So • Lno)^ — Sq . 

(41) 

We decouple Eq. ^T^i from Eq. then we use Eq. (|iT]) . 
Since we expect that 6r — 0(4), we find that at 2PN 
order 



Sr + Sr 



(So • Xf - (So • fi)^ 



(42) 



in agreement with Eq. (B13) of Ref. [45j. The solution 
of the homogeneous equation is simply 



<5r(t)hom = cos {Qt + (fir) , 



(43) 



where Ar and ipr are fixed by the initial conditions. 
Equation (^5)) describes possible oscillations due to the 
initial eccentricity of the orbit. This eccentricity occurs 
at the average orbital frequency and in principle can be 
removed as long as quasi-circular initial conditions are 
enforced. Note that Eq. is also consistent with the 
Newtonian result of Eq. (I3])I3 It is worth noting that 
Eq. ([3]) was derived as an expansion for a small eccen- 
tricity e, while in this section we never explicitly referred 
to e at all. As a matter of fact, we are dealing with quasi- 
circular orbits here, so that consistency between Eqs. ([3]) 
and (|43p should be expected. 

On the other hand, the spin-induced oscillations are 
described by the particular solution of Eq. (|^^ which 
reads 



<5rpart(i) 



1 



4A/2 



(So • A(f))2 - (So • fi(t))2 



(44) 



^ Note that Eq. l[3j docs not depend on ip because in Eq. ^ the 
radial separation at i = is picked to be equal to the semi-major 
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These osciUations are a signature of SS interactions since 
they depend on Sq, which enters i?ss- Once we know 
Sr — Srhom + ^fps^j-t, we solve Eq. (fSS)) for (5f2. Similarly to 
what we found above, the homogeneous solution accounts 
for the initial conditions, while the particular solution 
accounts for the oscillations induced by SS effects. 



n 



4M2 f2 



{So ■ xit)f - (So ■ mf 



(45) 



The above equation is also consistent with the Newtonian 
result of Eq. ([5]). 

So far, we have assumed the nonspinning dynamics to 
be Newtonian. If we included nonspinning PN correc- 
tions through 3PN order in the Hamiltonian, we would 
still find the particular solutions (|44| and p5)) . but in this 
case the oscillations would occur at a frequency which will 
differ from Eq. (j4ip because of nonspinning PN correc- 
tions. 

Using the previous results, it is straightforward to com- 
pute the time derivatives of the oscillations dr, Eq. (|33]), 
and 5Q, Eq. (gS]). They read 



6r{t) = Br sin (fit + f,.) 



5tl{t) = Bn sin {Clt + ipn 
02 



M2 f2 



(So-n(t)) (So-A(O), 
(So-n(t)) (So-A(t)). 



(46) 



(47) 



We note that when the spins are initially aligned or an- 
tialigned to Lno: the SS oscillations disappear, since 
in this situation So remains perpendicular to fi and A 
throughout the evolution. We sec that for both quanti- 
ties the time dependence of the SS term is 



(So-n(t)) (So-A(i)) = Csin(20t + 7), 



where 



C 



and 7 satisfies 



(So • fio)' + (So • Ao)=^ 



sin 7 ~ — sin(2a), 
cos 7 = — i cos(2a), 



(48) 



(49) 



(50) 



with cos a = So_l • no- Thus, the spin-induced oscillations 
occur at twice the average orbital frequency, and they 
can be neatly disentangled from the eccentricity induced 
by initial conditions which occurs at the average orbital 
frequency. 

Moreover, the amplitude of spin-induced oscillations is 
quite small. To place their amplitude into perspective, 
consider a binary with orbital eccentricity e. Taking a 
time-derivative of Eqs. ^ and ([S]), and comparing to 
Eqs. dUD and (gT]) we find Br = fCle and Bn = 2^2 e. 



Equating the amplitudes of the spin-induced oscillations 
with the amplitude of the eccentricity-induced term, we 
see that spin-induced oscillations dominate only for ec- 
centricities 



e < 



1 92 

2 M4 
152 



4 M4 



for Sr, 



for sn. 



(51) 



Numerical binary black-hole simulations typically start 
at a separation f/M 15, and in that case, spin- 
induced oscillations will dominate S-r and (5f2 only for 
e < 0.002S^^/M^ and for e < O.OOIS^ JM^, respec- 
tively. For maximally spinning black holes with least- 
favorable spin orientations, So±/M'^ = 1, so that even in 
this case spin-induced oscillations become relevant only 
at eccentricities of < 0.001. For smaller spins, their ef- 
fect is still smaller. We note that spin-induced oscilla- 
tions do affect SQ somewhat less than Sr, indicating that 
eccentricity-removal based on the orbital frequency will 
be preferable. 

Let us notice that were we including the precession 
of the spins, the characteristic frequency at which the 
spin-induced oscillations occur would change. This can 
easily be seen if we assume that the precession is mainly 
due to SO effects, with Si and S2 precessing about Ln 
at frequencies ili and ^2. In this case, using Eqs. (|15p . 
([TBI), we derive 



01 = 

02 = 



2/iO 

f 

2fiCl 



1 



3 7712 

4 mi 
3 mi 
4m2 



(52) 
(53) 



If in Eqs. (|48]) dH]) and (JSO]), we let So precess, we ob- 
tain that oscillations occur at frequencies given by linear 
combinations of H, Oi and ^2, namely 20 — Oi — O2, 
2(0 — Oi) and 2(0 — O2). For the binary black-hole 
evolutions considered in this paper, O1.2 ^ O, so the 
spin-induced oscillations occur at the frequency 20. 

Lastly, the results of this section could be extended to 
higher PN orders by including next-to-leading SO terms 
(2.5PN order ^MM) and SS terms (3PN order (H-iOl). 



D. Oscillations induced by leading SS terms: 
inspiraling dynamics 

In this section we compare the approximate analyti- 
cal predictions for Sr and (50 with the results obtained 
by numerically integrating the Hamilton equations (|13p 
and (|14p . including the radiation reaction force Ft. Since 
we actually want to extract the spin-induced oscillations, 
we need to remove the homogeneous part which is due 
to the eccentricity introduced by the initial conditions. 
In fact, in the presence of radiation-reaction the initial 
radial velocity has to be carefully chosen to guarantee 
that the binary moves along a quasi-adiabatic sequence 
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of spherical orbits, progressively shrinking. Those initial 
conditions have been worked out in the analytical PN 
dynamics at post-circular [6l] | and post-post-circular (62| 
orders. However, those initial conditions become more 
and more approximate if we start the evolution of the 
Hamilton equations at smaller and smaller separations. 
Moreover, from the study of the conservative dynamics in 
Sec. IH CI we understood that because of SS interactions 
it is impossible to have spherical orbits. To remove the 
oscillations at a frequency Cl (the homogeneous part), we 
perform a fit of the data with a function 

./oscin(^; Bat, Wfit, (ptit) = Bat sin(a;fitt + ^ptit); (54) 

where Wfit is close to O. We subtract the fitted function 
/osciii from the raw data sample, obtaining a residual that 
oscillates at a frequency 217, superimposed to the smooth 
numerical inspiral. The reason why we need to fit also 
the frequency uiat is that there is an ambiguity as to what 
value we use for O. In principle, this value should be the 
average orbital frequency, but a priori we can only use 
the initial value fto = il{0) because we do not have an 
analytic prediction for Q. This is also true for the value 
of f, which we replace with vq ~ r(0) throughout. We 
want to compare these residuals with analytical predic- 
tions based on Eqs. (|46)) and (|47)) . We use the expression 
of the Newtonian flux to derive the effect of the radiation 
reaction (RR) on the two quantities r and fl. Within the 
context of Newtonian dynamics, we have 



rRR 



4 256 , o 



-3/4 



where rg = r(0) and 



1/4 



(55) 



(56) 



A similar expression can be found also for JIrr . Consid- 
ering that in a quasi-circular inspiral 



M, we find 



2n 



in 



RR 



RR 



3 , ,1/2 rRRjt) 
5/2 ^ 



(57) 



(58) 



'rr(0 

Therefore our analytical predictions will be given by 

r-pi-cdit) = r-RRit) + (5fpart(i) (59) 
flpred {t) = flRR {t) + <5f2part (0 , (60) 



where 5rpart(^) and (5f2part are given by the second term 
in the RHS of Eqs. gg and (gZl). 

Figure [5] shows the results for a particular black-hole 
binary with mass ratio q = m^/mi ~ 2 and maximal 
spin magnitudes jSil/rnf = |S2|/m2 ~ 1. In these plots 
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FIG. 2. Spin induced oscillations in PN model. We com- 
pare two PN calculations of spin-induced oscillations. The 
dashed lines are the predictions from Eqs. 1)59^ and H60|) . with 
the two hues differing by whether Ln is held constant ( "non- 
precessing") or evolving ("precessing"). The sohd line repre- 
sents a solution of the full PN equations of motions, with 
the homogenous oscillations fitted and subtracted. (Here, 
mass-ratio q — 2, with maximal spins of initial orientations 
(^1 = 7r/3,(^i = 0) and {62 ~ 27r/3, (^2 ~ 7i"/3) at the initial 
orbital frequency MHq = 0.015, that is an average initial or- 
bital period P = 418M, and the dotted curve represents the 
solutions of the PN model including inspiral motion). 



the initial orbital frequency is floM — 0.015 which corre- 
sponds to a period P = 418M. As can be seen, the raw 
data (dotted lines) are dominated by the eccentricity at 
the orbital frequency, while the residuals (solid lines) os- 
cillate at twice the orbital frequency, as expected. As 
far as the amplitude of the residual oscillations is con- 
cerned, we see that it is compatible with that of the pre- 
dictions computed using Eqs. ([SH]) . even though the 
agreement is not striking (red dashed lines, called non- 
precessing). In fact, we find that the effectiveness of the 
removal procedure of the homogeneous part is deeply af- 
fected by the value of Wfit ■ Numerical studies have shown 
that differences of a few percent on ujat can completely 
alter the residuals. Tweaking by hand the value of ujat 
instead of using the best fit value can lead to a much 
better agreement on the amplitudes, at least for the first 
cycle. Note that in Fig. [5] no such ad hoc tweaking is 
used. Also, the fact that the predictions quickly get out 
of phase with respect to the residuals is mainly due to 
the assumption made in Sec. Ill CI of keeping the spins 
constant or non-precessing, that is using the evolution 
of the triad {n. A, Ln} specified by Eqs. p3|) - ((35|) . By 
contrast better phase agreement can be obtained numer- 
ically if we use the time-evolution of the spins and of 
the reference triad (blue dashed lines, called precessing) 
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or analytically if we had considered a reference triad in 
which the precession of the orbital plane and spins were 
taken into account. 

Another interesting feature is the relative importance 
of the spin-induced oscillations with respect to the ec- 
centricity induced by the initial conditions. Both types 
of oscillations are showing up in the raw data of f, while 
in the case of ft we only see the oscillations due to the 
initial conditions. This can be explained by our analyt- 
ical predictions. Using Eqs. ([2]), ([5]), we write for the 
eccentricity 



.,NS 



or 



„NS 



2n^ ' 



(61) 



where Br and arc the amplitudes of the oscillations 
of the homogeneous solutions. We want to keep distinct 
notations for f and O, even though at Newtonian level 
= e^s ^ ^ therefore \Br\ = f\Bn\/2n . If we 
now call Cr and Cn the SS amplitudes, namely 



Cr 



nc 



Cn = - 



M2 f2 ■ 



(62) 



we parametrize the spin-induced oscillations in terms of 
a spin-induced "eccentricity" , 



lai \c\ 



OS _ ICnl _ \C\ 



^ 2n^ 2M^ f2 ■ 

(63) 

We have that the relative ratio is 



oSS 



1 



(64) 



so it is expected that the significance of the spin-induced 
oscillations is smaller for Cl. 



III. ITERATIVE ECCENTRICITY REMOVAL 
IN PRESENCE OF SPINS 

In the preceding sections, we showed that the PN 
Hamiltonian with leading SS terms predicts oscillations 
with two distinct periods: the orbital period, with am- 
plitude and phase depending on initial conditions that 
can genuinely be associated with orbital eccentricity; and 
half the orbital period, independent of the initial con- 
ditions and caused by spin-spin couplings. We further- 
more showed that these spin-induced oscillations are sup- 
pressed in (l as compared to r [see, e.g., Eq. (|64| ]. 

Our task is to find initial conditions that remove or at 
least minimize the oscillations caused by eccentricity. As 
in earlier work, we shall begin with some trial initial con- 
ditions, evolve the binary for about two orbits, analyze 
the motion of the black holes, and then correct the initial 
conditions. To exploit this suppression of spin-induced 
oscillations in f2, we will derive updating formulae based 
on f2. 



A. Updating formulae 

The basis for the updating formulae arc the Newto- 
nian expressions for distance r and orbital frequency Q, 
(Eqs. © and ©) 



rN{t) = r [1 - e sin0(<)] , 
VtN{t)^Vt [l + 2e sin0(i)] 



(65) 
(66) 



where 0(i) is the phase of the radial oscillations. General 
relativistic periastron advance will cause (j){t) to deviate 
from the orbital phase. Taking a time-derivative, we find 



f^V ~ e LO COs{LUt + 0o)j 

r^AT = 2Q e Lo cos(ajt + (po), 



(67) 
(68) 



with uj = {d(j)/dt){0), 00 = ^(0). 

Let us now consider a compact binary inspiral starting 
at t = at initial separation vq, with orbital frequency 
Hq and radial velocity t-q. We take r(t) or (l{t) from 
a general relativistic inspiral, and fit it with functional 
forms 

rNR{t) = Sr{t) + Br COs{uJrt + (j)r) , (69) 

nM-t) = Sn{t) + Bn cos{Lunt + (70) 

The subscripts r and 51 indicate whether the fit was per- 
formed on tnr or f2NR, and we will use a bullet • in the 
subscript to represent either r or fl. The first part of 
each fit, iS,, is a non-oscillatory function that captures 
the radiation-reaction driven inspiral, whereas the oscil- 
latory piece captures the orbital eccentricity. We neglect 
spin-induced oscillations. The precise functional form of 
S, is important, and sometimes it is advisable to include 
a quadratic term Ct^ within the argument of the cosine. 
We comment on these considerations below in Sec. IIIIBI 
Equation shows that at i = 0, orbital eccentricity 
contributes Tocc.o = Br cos 0^ and rccc,o = —BrUJr sin (j)r 
to the radial velocity and acceleration. Our goal is to 
now modify the initial data parameters fo and Qq such 
that focc,o and fecc.o vanish. The radial velocity is a free 
parameter of the initial data, so fo.now = 'f'o + Af, where 



Af = -r-ccc.O = -i?rCOS0r- 



(71) 



To utilize our information about the radial acceleration 
^ecc.o we recall that for the Newtonian Hamiltonian we 
have 



Pr 1 dH^ 



/.J /.J dr fi 



L2 M ^, M 
— rff 



(72) 



A small change Oq — > ^o.new = ^o + Ail therefore changes 
the radial acceleration by Af = 2roiloAil. This change 
cancels f ccc o when 



All 



2roilo 



(73) 



9 



■l-H 




•a 



— Iter 0, e=0.00328 

— Iter 1, e=0.00093 

— Iter 2, e=0.00035 




FIG. 3. Eccentricity removal based on Q applied to a PN 
model. Shown is the initial orbital evolution, and two iter- 
ations of eccentricity removal. Parameters of the black-hole 
binary: mass-ratio q = 1, spins of magnitude Xi = X2 = 0.8 
with initial orientations {8i = 0,(l>i ~ 0) and {02 — (j)2 ~ 
0) with respect to Lno and initial AIQq = 0.0315, that is an 
initial orbital period P = 200M. 



Equations ([7T|) and ((75|) are one version of the updating 
formulae for the eccentricity removal based on the sep- 
aration coordinate. A sometimes more effective formula 
is presented below in Eq. (|77p which in earlier numeri- 
cal work [Hill [Hi was applied to the proper separation 
between the horizons of the black holes. 

A convenient way to derive updating formulae based on 
Cl{t) begins by noting that the ratio of the amplitudes of 
oscillations in Eqs. ^ and ^ is -f/(20). Therefore, 
we obtain the desired updating formulas by replacing B^. 
with -fBn/{2fl): 



Ar 



2rJo 



cos 0n, 



iXil = —TT- smc 



(74) 
(75) 



[A sometimes more effective replacement for Eq. ([75]) is 
presented below in Eq. ((78)) .] 

In Fig. [3l we present three steps of this eccentricity re- 
moval procedure. We use the PN-expanded Hamiltonian 
with non-spinning terms up to 3PN order [sol |63| , and 
spinning terms up to 2PN order [i^ [s^j . The radiation- 
reaction effects arc included through 2PN order as in 
Ref. 

We indicate in the legends the value of the eccentricity 
estimated from the fitted amplitude of the oscillations, 
once the smooth inspiral has been removed. Note that 
the initial orbital period is about 200M. At the O"' step 
the plots are showing the evolution of the binary system 



with initial conditions determined according to the pro- 
cedure outlined in Ref. [6l| , leading to the presence of an 
initial eccentricity. From the plots, we clearly see that we 
go from a situation dominated by the homogeneous oscil- 
lations occurring at the average orbital frequency (step 
0) to the situation in which only spin-induced oscillations 
occurring at twice the average orbital frequency are vis- 
ible (step 2). 

The configuration considered in Fig. [3] is close to 
merger, where the rapid evolution of the orbit makes it 
more difficult to apply eccentricity removal. In the next 
section we discuss how to improve the convergence of the 
iterative procedure. 



B. Practical considerations 

Unfortunately, iterative eccentricity removal is sensi- 
tive to a variety of effects which are not immediately 
obvious. Without sufficient care, iterative eccentricity 
removal converges slowly, or not at all. In this section, 
we describe important details for the effective and prac- 
tical application of the eccentricity removal, as well as 
diagnostics that allow users to evaluate whether the ec- 
centricity removal proceeds optimally. 

The fits in Eqs. ([BQ]) and ([70)) are used to compute the 
values of the oscillating part B, cos(w,t -|- (/>,) at t = 0. 
Therefore it is crucial that the function S, that is in- 
tended to fit the inspiral portion docs not fit this oscil- 
latory piece. Initially, we used a polynomial for S*,, but 
sometimes, especially for shorter fitting intervals, such 
a polynomial picks up a contribution of the oscillatory 
piece and results in an unusable fit. Therefore we have 
constructed more robust fitting functions that cannot 
capture oscillations. Our current preferred choice is 



k-l 



ll/8-n/4 



(76) 



n=0 



with free parameters Ak and Tc- The functional form 
and the exponents are motivated by PN inspirals, and 
we keep either fc = 1 or fc = 2 terms of this expansion 
(for fc = 2, we use the same Tc in both terms). 

Another crucial ingredient for a reliable fit is a suit- 
ably chosen fitting interval. This interval needs to cover 
enough oscillations to break degeneracies among the fit- 
ting parameters. However, if it becomes too long, the 
evolution in the inspiral part will be harder to capture 
with St and the quality of the fit will deteriorate. 

Finally, the fit is used to compute Tecc.o and recc.Oj 
which are quantities at i = 0. It is desirable that the 
fitting interval starts as close to < = as possible to 
minimize extrapolation from the fitting interval to t = 
0. However, a numerical evolution relaxes in its early 
stages due a quasi steady-state, and features during this 
relaxation need to be excluded from the fitting interval. 

A good means to ensure a satisfactory fit is to perform 
several fits, and ensure that the results are consistent. 
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FIG. 4. Visualization of the eccentricity removal performed 
in Fig. [3] in the Qq-clo plane. This plot summarizes a large 
amount of diagnostic information which can be utilized to 
ensure reliability of the eccentricity removal procedure (see 
main text). The symbols B, and uj have a subscript Q 
suppressed for clarity. 



We perform four distinct fits to 17, where we change the 
order fc = 1,2 of the inspiral component, Eq. ([76]) . and 
where we change the order of the polynomial within the 
cosine in Eq. ([70]) between I = 1 as shown in Eq. ((70)) and 
I ~ 2 (i.e. adding a quadratic term C,t^). In addition, 
we vary the location and length of the fitting interval and 
check that the updates Sr and Ail are unaffected. 

Figure |4] demonstrates a useful way to visualize and 
assess the effectiveness of iterative eccentricity removal. 
This plot shows the plane of the initial-data parameters 
^0,00, with ciQ = ro/ro- The solid symbols correspond 
to the three PN evolutions shown in Fig. [31 The lines 
emanating from each of the symbols denote different fits 
based on this particular evolution. The different fits are 
denoted Fk cos;, where fc = 1, 2 denotes the fitting order 
in 5", and / denotes the order of the polynomial inside the 
cosine. Each of these lines ends at the predicted improved 
parameters rio,new, ao,new Clustering of these lines, and 
convergence of the symbols indicates that eccentricity re- 
moval is proceeding well. As can be seen by comparing 
the solid dark green and dashed red curves, the order k 
of the fitting function for the smooth part Sq has almost 
no impact on the updated parameters llo,ncw7 ao,ncw h^ 
this case. However, using a quadratic polynomial inside 
the cosine (l = 2) significantly improves the quality of 
the f2o~update. 

FigurelUcan also be used to assess the potential quality 
for different updating formulae. While we have kept the 
orbital frequency VIq separate from the eccentricity oscil- 
lation frequency ujq, for Newtonian orbits both frequen- 
cies agree. Therefore, our Newtonian motivation does 
not provide a means to choose whether to include extra 



powers of flo/ujQ. Specifically, we could replace Eqs. ((73)) 
and ()75p by either 

AO = sm(t)r, (77) 

2ro 

or 

Af^ = -^sin0n, (78) 

for updates based on r{t) and (i{t), respectively. The 
predictions of the updating formula Eq. (|77p are shown in 
Fig. ID as filled blue triangles. It is obvious that Eq. ((77)) 
predicts an updated Oq significantly closer to the final 
best value for Qq, even when applied to Iter 0. Therefore, 
to summarize the discussion in the preceding paragraphs, 
for most effective eccentricity removal we recommend the 
fit of the form F2COS2 combined with Eq. ((77)) . 

Finally, we discuss several diagnostics that can help to 
assess the quality of eccentricity removal, and which are 
included next to each symbol in Fig. H) The first diag- 
nostic is the estimated eccentricity e, which should be 
monotonically decaying. The second diagnostic is sin (pQ . 
As can be seen from Eqs. ((7^ and ((75)) . the angle 0n 
parametrizes the relative importance of the flo and clq up- 
dates. For sin ~ 1 , the whole weight is carried by the 
fto update. This is the case here, and indicates that the 
starting value for do was already very good, and that the 
apparent inconsistent predictions for do, new will not have 
an adverse impact on the eccentricity fitting procedure 
(note that all fits predict consistent values for rio,ncw)- 
The third diagnostic is the ratio of the root-mean-square 
residual of the fit, res, to the amplitude of the oscilla- 
tory part, Bq. When rcs/B <C 1, then (l has indeed 
the assumed form Eq. ()70p . a prerequisite for eccentric- 
ity removal. When ics/B ^ 1, we can no longer isolate 
the oscillatory piece, and eccentricity removal ceases to 
be effective. The final diagnostic is the ratio of frequen- 
cies of radial oscillations wq to orbital frequency Qq. For 
a good fit. wo/Oo should be somewhat less than unity, 
where the deviation from unity is caused by pcriastron 
advance. For moderately small eccentricities, this ratio 
should furthermore be independent of the precise value 
of eccentricity. This is indeed the case for "Iter 0" and 
"Iter 1", but "Iter 2" results in a questionably small ratio, 
which furthermore differs from the values for iterations 
and 1. Again, an indication that we cannot further 
proceed with eccentricity removal. 

IV. APPLICATION TO FULLY NUMERICAL 
BINARY BLACK-HOLE SIMULATIONS 

We now apply the method outlined in Sec. lIIII to reduce 
the initial eccentricity of single-spin and double-spin pre- 
cessing binary black-hole simulations. We compare the 
periodicity in the oscillations of the orbital frequency and 
the proper horizon separation to the PN predictions de- 
scribed in Sec. |TT] and also to the periastron-advance re- 
sults of Ref . [J^l . Finally, for one binary configuration, we 
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also extract the I = 2, m = 2 mode of the gravitational 
waveform and investigate the presence of spin-induced 
oscillations in its phase and frequency. 



A. Numerical methods 

Binary black hole initial data is constructed using the 
conformal thin sandwich formalism J64, J65| and quasi- 
equilibrium boundary conditions (28l . l66l. l67ll . incorporat- 
ing radial velocity as described in Ref. [3j| . The resulting 
set of five nonlinear coupled elliptic equations is solved 
with multi-domain pseudo-spectral techniques described 
in Ref. [i^ . As in earlier work, we choose conformal flat- 
ness and maximal slicing. To obtain desired masses and 
spins, we utilize a root-finding procedure to adjust freely 
specifiable parameters in the initial data j40| . 

Thus, a binary black-hole initial data set is determined 
by the mass-ratio, the spins of both black holes, and co- 
ordinate separation d between the coordinate centers of 
the black holes, orbital frequency iloi and radial velocity 
fg = aod^ where dg is the dimcnsionless expansion factor. 
For fixed d, eccentricity removal consists of finding values 
for flo and do that result in sufficiently small eccentricity. 

The constructed initial data are evolved with the Spec- 
tral Einstein Code SpEC [69l |. This code evolves a first- 
order representation [70| of the generalized harmonic sys- 
tem j7ll473{ and includes terms that damp away small 
constraint violations The computational do- 

main extends from excision boundaries located just inside 
each apparent horizon to some large radius. No bound- 
ary conditions are needed or imposed at the excision 
boundaries, because all characteristic fields of the system 
are outgoing (into the black hole) there. The boundary 
conditions on the outer boundary [t^, [iE IZ^ are de- 
signed to prevent the infiux of unphysical constraint vio- 
lations 
diation 



77 

84 



and undesired incoming gravitational ra- 
85| , while allowing the outgoing gravitational 
radiation to pass freely through the boundary. Interdo- 
main boundary conditions are enforced with a penalty 
method [Hiai. 



B. Eccentricity removal based on orbital 
frequency: single-spin binary black hole 

In this section we re-visit eccentricity removal for the 
configuration considered in the introduction and Fig. [T] 
The binary has a mass-ratio of toi/to2 = 1.5, and only 
the larger black hole carries spin, namely xi = 0-5 with 
initial spin direction in the orbital plane pointing exactly 
away from the smaller black hole. Note that spins in 
the orbital plane maximize spin-induced oscillations [see, 
e.g., Eqs. (HH) and (US])]. The initial coordinate sepa- 
ration between the holes is d = 16M, orbital frequency 
MQq = 0.0142, and do = — 5 x 10~^. These parameters 
were determined from the so-called TaylorTS FN approx- 
imant for non-spinning binaries [l6j . 



30 



20 



.a 



10 



-10 



0.8 







- IterO: 


e~0.04 






— Iterl: 


e-0.01 






OG Iter 2: 


e~0.0014 






>^ Iter 3: 


e-0.0001 






Iter 4: 


e-0.00006 





1 

t/p 



FIG. 5. Eccentricity removal based on time derivative 
of the orbital frequency dQ,/dt, applied to a single-spin 
precessing binary black hole with the same initial parameters 
as in Fig. \T\ Shown is Cl vs. time in units of initial orbital 
period P = 442M for the initial run (based on FN param- 
eters) and four eccentricity-removal iterations. The ampli- 
tude of spin-induced oscillations is several orders of magnitude 
smaller than in Fig. [T] and becomes only visible in Iter 4. 
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FIG. 6. Convergence of the eccentricity-removal pro- 
cedures in the (f2o5<io) plane. Blue circles: eccentricity 
removal sequence of Fig. (5] Red squares: eccentricity removal 
sequence shown in Fig. [T] The number next to each symbol 
gives the eccentricity of the respective evolution. The inset 
shows an enlargement of the boxed area. The eccentricity- 
removal procedure based on the orbital frequency keeps con- 
verging until the fourth iteration, while the one based on the 
proper separation fails to converge any further beyond the 
second iteration. 
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FIG. 7. Radial velocity between the black holes for the same series of evolutions shown in Fig. (5] Left: Derivative of proper 
separation between the apparent horizons s. Right: Derivative of coordinate distance between centers of the apparent horizons 
f. The time-axis is given in units of the initial orbital period P — 442Af. Proper separation s exhibits large spin-induced 
oscillations, whereas r shows spin-induced oscillations of similarly small amplitude as in Fig. [S] Note that the inset in the right 
panel is ten times more magnified than in the left panel. 



Orbital frequency, both in the initial data and in the 
subsequent evolution, is defined by the coordinate mo- 
tion of the center of the apparent horizons. Let Ci(t) be 
the coordinates of the center of each black hole, and de- 
fine their relative separation r{t) = Ci{t) — C2{t). The 
instantaneous orbital frequency is then computed as 



and ft is defined as the magnitude of ft. All these cal- 
culations are performed using standard Euclidean vector 
calculus. 

We start the first run using PN initial conditions for 
the orbital frequency and radial velocity and evolve the 
binary black hole for about two orbits. From the or- 
bital frequency we measure an eccentricity e ~' 0.04, and 
Eqs. ([71]) and ([75]) give improved values for Hq and ctQ. 
Evolution of the initial data computed from these im- 
proved values is labeled "Iter 1" in Fig. [5l and reduces 
the eccentricity to about 0.01. The same procedure is 
then repeated three more times. For Iter to Iter 2, we 
exclude t < lOOAf from the fit. For Iter 3 the variations 
in fl are so small that numerical noise is dominant for 
about half an orbital period, and we exclude t < 250M 
from the fit. The final eccentricity in the orbital fre- 
quency is e ^ 6 X 10~^. 

In Fig. |6l we show how the initial orbital frequency flo 
and the radial expansion factor do converge to the final 
values (minimal eccentricity). The blue circles indicate 
the successive iterations for the successful eccentricity re- 
moval based on fl (see Fig. [5]). Note that the parameters 
{ilo,ao) converge well for all iterations. In contrast, the 
red squares denote the unsuccessful eccentricity removal 
based on proper separation s (see Fig. [Ij . Starting with 



the third iteration, the updated values of the orbital fre- 
quency and expansion radial coefficient move away from 
the line of minimum eccentricity, with an increase of ec- 
centricity from 0.001 to 0.003. All eccentricity estimates 
shown in this figure are computed from fi, even when 
eccentricity removal is based on s. This allows us to 
measure eccentricities e < 0.01, which would not be pos- 
sible from s, because the latter is dominated by large 
spin-induced oscillations. 

The absence of spin-induced oscillations in Fig. [S] is 
striking, especially when compared to Fig. [1] Spin- 
induced oscillations are visible in Fig. [S] only at e < 
lO""'. For the runs with larger eccentricity (Iter 0-3), 
eccentricity-induced oscillations dominate with a period 
somewhat larger than P (somewhat larger because of 
periastron-advance [13 )• 

We shall now investigate spin-induced oscillations in 
the numerical-relativity simulations in more detail. First, 
by comparing the time-derivatives of orbital frequency 
n, proper separation between horizons s, and coordinate 
separation between centers of apparent horizons r. Sub- 
sequently, by comparing their amplitude and frequency 
with PN predictions from Sec. |IT1 

Figure [7] shows time-derivatives of proper separation s 
and coordinate separation r for the evolutions shown in 
Fig.[SJ Spin- induced oscillations are already noticeable in 
s for Iter 1 with e = 0.01. These oscillations dominate for 
Iter 2-4, i.e. e < 0.0014. In contrast, the coordinate dis- 
tance f is less susceptible to spin-induced oscillations. In 
the right plot of Fig. [71 spin-induced oscillations become 
apparent only for eccentricities of 10~* or smaller. The 
spin-induced oscillations in r are smaller by a factor of 
almost 20 than those in s. 

When comparing Iter 3 and Iter 4 between Fig. [5] and 
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FIG. 8. Spin-induced oscillations of the lowest-eccentricity 
numerical simulation (Iter 4 of Fig. [5]) in comparison with the 
PN calculations of Sec. HTCl 



the right panel of Fig. [71 one notices that (l shows slightly 
less pronounced spin-induced oscillations. That is consis- 
tent with the PN calculations, where Eq. (|64|) predicted 
that spin-induced oscillations in f2 should be suppressed 
by a factor of 2 relative to those in f. When comparing 
Iter 4 between Fig. [5] and the right plot of Fig. [71 we 
find that the spin-induced oscillations in O and r are in 
phase. This is again consistent with the PN prediction, 
where the last terms of Eqs. (|46p and (|47| have the same 
sign. The phase of spin-induced oscillations in f2 and r 
differs from the effect of orbital eccentricity: for an ec- 
centric orbit, the orbital frequency is maximal when the 
separation is minimal, and therefore and f arc out of 
phase (see Iter and Iter 1 of Figs. [Hand [7]). 

We have shown in Sec. Ill CI that the PN Hamiltonian 
predicts spin-induced oscillations: Equations (|^5|) - P^ 
contain an oscillatory component at twice the orbital fre- 
quency with amplitudes 



Asr 



q2 c2 



(80) 
(81) 



Figure [HI shows numerical data for the lowest-eccentricity 
numerical simulation (Iter 4 from Fig. [SJ. These numer- 
ical data are compared with the prediction of the PN 
equations. The PN calculation reproduces very accu- 
rately the amplitude of spin-induced oscillations in the 
numerical- relativity simulation for Q and r. By contrast, 
the spin-induced oscillations in s are larger by a factor 
~ 20 than those in r. This can be due to deformations 
of the apparent horizons due to spin effects. Finally, we 
notice that a small amplitude oscillation of the numerical 
data on the orbital time-scale is also visible, correspond- 



ing to the small, but non-zero eccentricity e = 6 x 10 ^ 
of the numerical simulation. 



Oscillations in the (2,2) mode of the 
gravitational wave 



In Sec. IIVBI we found spin-induced oscillations in the 
coordinate distance of the black holes and the orbital fre- 
quency, consistent with PN predictions. We now inves- 
tigate the gravitational radiation emitted by this binary. 
Specifically, we extract the I = 2,m = 2 mode of the 
gravitational waveform in the inertial frame where the 
binaries are initially placed along the a;-axis and the ini- 
tial angular momentum is parallel to the z-axis. We com- 
pute phase and frequency for the waveforms extracted at 
extraction radii R = 13GM and R = 220M. 

Spin-induced oscillations represent a physical effect 
independent of orbital eccentricity. Nevertheless, the 
concept of eccentricity estimators [13] will be very use- 
ful when discussing spin-induced oscillations, because 
it removes overall secular trends (especially in the 
gravitational- wave phase), and because it makes it easy 
to relate the amplitude of oscillations to an "equivalent 
eccentricity." As Ref. [i^l , we define e^^^ (t) 



it) 



ffit 



it) 



(82) 



where 0nr(O is the gravitational-wave phase of the (2,2) 
mode and (/)fit(0 is the quasi-circular polynomial fit of 
the gravitational-wave phase [see Ref. [421 for more de- 
tails]. Similarly, using the gravitational- wave frequency 
of the ( 2,2) mode and its quasi-circular polynomial fit as 
in Ref. [42| , we define the eccentricity estimator Ci^^w (0 



At) 



jONRjt) - ujmjt) 
2u;fit(<) 



(83) 



In Fig. [9l we plot the eccentricity estimators e^^^ {t) 
(upper panel) and e^^,^ (t) (lower panel) for the two ex- 
traction radii R = 130M and R — 220M versus the 
retarded time t — R*, where R* is the tortoise-coordinate 
radius defined as 



i?*.i? + 2Mln(A 



1 



(84) 



where Af = 1 is the total mass of the initial data. Quite 
interestingly, the plots show oscillations happening at 
twice the orbital frequency. The magnitude of the os- 
cillations in e0Q„(t) or ei^^^ (t) is ~ 10"^, ahhough the 
eccentricity in the orbital frequency has been reduced to 
^ 6 X 10~^ (see Fig. [5]). We note that the amplitude 
of the oscillations at twice the orbital frequency does not 
depend on the extraction radius, suggesting that they are 
gauge invariant. 

We also compare these numerical result with what is 
predicted by the PN model. For the orbital evolution 
we use the model Hamiltonian ((51), where SO and SS 
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FIG. 9. Eccentricity estimator for the gravitational- 
wave phase and frequency for the final eccentricity- 
removal iteration of Fig. (5] The upper panel shows the 
gravitational-wave-phase eccentricity estimator for two ex- 
traction radii versus the retarded time, the lower panel the 
eccentricity estimator computed from the gravitational wave 
frequency. In both panels we clearly see oscillations at twice 
the orbital frequency. We also show the PN result for the 
restricted waveform. 
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FIG. 10. Eccentricity removal for different spin- 
configurations. We illustrate of how the eccentricity is re- 
duced to very low values when the iterations are applied to 
the orbital frequency. Shown are three configurations with 
S-l/M? = 0.5, S2 = and different spin directions 61, and one 
configuration with Si /Mi =82/ M2 =0.5, with initially two 
orthogonal spins both tangent to the orbital plane. For all 
cases, the mass-ratio is m\/m2 = 1.5. The run shown in green 
triangles was discussed in detail in Figs. [5] to [O] 



interactions are included through 2PN order, nonspin- 
ning effects through 3PN order, and radiation-reaction 
effects through 2PN order. As to the analytical model, 
we employ the waveform derived in Ref. [88|, where the 
precession of the orbital plane and the spins of the black 
holes are taken into account through 1.5PN order. In 
particular, we compute the estimators e^j,^ and e^^^^ 
using Eqs. (4.15), (4.16a) in Ref. [Ill for the /122 mode 
with the amplitude computed at lowest order in v /c (re- 
stricted waveformj^. This means that the precession of 
the orbital plane is considered only in the gravitational- 
wave phase, but not in the amplitude. 

Before computing the PN eccentricity estimators, we 
perform an alignment between the phase of the PN /122 
and the phase of the numerical-relativity '3/4. To do that, 
we follow the procedure outlined in Sec. IIIA of Ref. @. 
This alignment is obtained over a time window of lOOOM 
(in the range 500 A/ < t < 1500M), and it amounts to 
a time-shift and a global offset in the PN phase. The 
result is shown as solid black curves in Fig. [9l We see 
a qualitative agreement between numerical-relativity re- 
sults and the restricted PN model for the oscillations at 



* Since the numerical-relativity (2,2) mode is computed using as 
2-axis the direction perpendicular to the orbital plane, we apply 
the Wigner rotation to the restricted h22 of Ref. ISSlI and keep 
only the lowest-order term in v/c. 



twice the orbital frequency in e^^^ and e^^^. However, 
the numerical-relativity e,^^^ also shows oscillations at 
the orbital frequency which are absent in the restricted 
PN waveform. We find that oscillations at the orbital 
frequency can be generated in the PN model of Ref. [Ill 
if we included higher order PN corrections in the ampli- 
tude of the (2,2) mode [see Eq. (4.16a) in Ref. [H]. Such 
oscillations cannot be iterated away by our procedure, 
even in principle, since the removal algorithm concerns 
the orbital dynamics and they rather appear as a physical 
effect of the waveform. The upper panel of Fig. [H] shows 
comparatively large oscillations at period ^ P; because 
the PN model predicts modes at this frequency, these 
oscillations cannot be used to compute orbital eccentric- 
ity. A further analysis of the inclusion of higher-order 
PN corrections is warranted. Wc prefer to postpone such 
an analysis to be able to test against a larger sample of 
numerical-relativity waveforms. 



D. Eccentricity-removal for generic binary black 
holes 

In the previous sections we studied our new 
eccentricity-removal procedure in detail for one test-case: 
a binary with only one non-zero spin, and with mass- 
ratio mi/m2 = 1.5. We now test the procedure for 
other binary configurations with the same mass ratio. 
We consider two further configurations where only the 
large black hole carries spin, parametrized by the angle 
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01 between the orbital angular momentum and the spin 
axis of the first black hole. In the previous sections, we 
considered 6i = tt/2, and now we extend to 9i = 0, 7r/6, 
and 01 = 0. The first of these new cases is non-precessing 
and verifies that eccentricity removal based on D, works 
effectively for non-precessing systems. We also consider 
a binary where both black holes carry spin, with initial 
spin-direction in the orbital plane {9i = 02 = 7r/2), Si 
parallel to the initial separation vector between the black 
holes, and S2 normal to the separation vector. (All spin- 
ning black holes have dimensionless spin-magnitude of 
0.5.) Figure [TUl demonstrates the effectiveness of the ec- 
centricity removal procedure based on and Eqs. (j74p 
and ([75|) . In all cases, the eccentricity is reduced to less 
than 10~^ in four iterations. 

The number of required eccentricity removal iterations 
depends on the quality of the guess for fio and do for the 
first iteration. Once eccentricity removal has been per- 
formed for several different configurations, wc expect to 
be able to interpolate between configurations, to improve 
the initial guess substantially. 



V. CONCLUSIONS 

The removal of the initial spurious orbital eccentricity 
in binary black-hole simulations is quite challenging, and 
it becomes more complicated in the presence of spins. As 
predicted by PN theory, and worked out in Sec. |IT1 spin- 
spin interactions (notably Si Si , S2 S2 and Si S2 terms) 
and precession induce oscillations in the binary radial 
separation and orbital frequency. These spin-induced os- 
cillations are also present in the gravitational radiation 
emitted by the binary, and their frequency is close to 
twice the average orbital frequency. In Sec. lIVI we confirm 
the presence of spin-induced oscillations in fully numer- 
ical simulations of black hole binaries. The agreement 
between the numerical simulations and PN calculations 
is as good as can be expected given the low order of the 
PN calculations, and the differing coordinate gauges. 

Spin-induced oscillations can be distinguished from os- 
cillations caused by orbital eccentricity by their charac- 
teristic frequencies. Moreover, the amplitude of spin- 
induced oscillations is quite small, so that it becomes 
visible only at small eccentricities, as can be seen from 
Eq. (|51|) : At separations relevant for numerical simula- 
tions, spin-induced oscillations dominate orbital eccen- 
tricity only for e < 0.001, even for maximal spins in 
the least favorable orientation (parallel spins in the or- 
bital plane). The amplitude of spin- induced oscillations 
is proportional to Sq j^, so that for spin S/A'P ^ 0.5 as 
considered here, spin-induced oscillations become visible 
at orbital eccentricity e ~ 10~^. 

Spin-induced oscillations affect the orbital frequency 
derivative (l less than the radial velocity f. Therefore, 
we develop in Sec. IIIII a new eccentricity-removal pro- 
cedure based on the derivative of the orbital frequency, 
and apply it to PN inspirals. Subsequently, we success- 



fully apply the eccentricity removal procedure to fully 
numerical binary black hole evolutions to achieve eccen- 
tricities smaller than 10""'. At this residual eccentricity, 
spin-induced oscillations begin to dominate over orbital 
eccentricity oscillations, as expected from our PN cal- 
culations. In Sec. IIV D[ we tested the new eccentricity- 
removal procedure on fully numerical binary black hole 
simulations for several different spin configurations. 

The computational cost for eccentricity reduction de- 
pends on the number of eccentricity removal iterations. 
Great care is necessary when performing the fitting, in 
order to achieve a large reduction in eccentricity per it- 
eration. Section IIII Bl gives guidance to improve these 
fits. With a good initial guess of flo, do for the 0-th iter- 
ation, one can start eccentricity removal from an already 
small initial eccentricity. As the number of data-points 
increases, we expect to be able to compute a better initial 
guess by interpolating between already computed low- 
eccentricity binary black- hole configurations. 

Perhaps surprising, our present study indicates that ec- 
centricity removal should not be based on the proper sep- 
aration between the apparent horizons. These new find- 
ings supersede the practice of earlier papers [l^ 
to base eccentricity removal on proper separation rather 
than coordinate separation to take advantage of reduced 
numerical noise. As apparent in Figs. [7] and [8l spin- 
induced oscillations in s are about 15 times larger than 
in r. Therefore, eccentricity-removal based on s will fail 
at 15 times larger eccentricity than using f, and at 
^ 30 times larger eccentricity than for Cl. A likely cause 
for the unsatisfactory behavior of s lies in the deforma- 
tion of the apparent horizons due to spin. For spins with 
a component in the orbital plane, the "bulge" of the ap- 
parent horizon rotates through the line connecting the 
two black holes as the black holes orbit each other. Ear- 
lier work [1^, [2^ [11] considered spins aligned with the 
orbital angular momentum, where this effect is absent; 
in those cases use of s was in order — but for precessing 
binaries, use of s is not advisable. 

Even when the orbital frequency indicates e < 10~^ 
for a fully numerical binary black-hole simulation, the 
extracted (2.2) mode of the gravitational radiation still 
shows oscillations at the orbital frequency in the wave 
phase with amplitude ~ 10"'^. While future work is nec- 
essary for a detailed understanding, the PN model pre- 
dicts oscillations in the GW at the orbital frequency, and 
therefore, one cannot use the gravitational waveforms 
to estimate orbital eccentricity for precessing binaries. 
The wave phase and frequency of the NR simulation also 
shows oscillations at twice the orbital frequency with am- 
plitude ^ 3 X 10^* which arc qualitatively reproduced 
by the restricted PN model of Rcf. Wc postpone 

the study of the details of these features in the gravita- 
tional waveform to future work. Quite interestingly, it 
proves, for this particular binary configuration in which 
only one hole spins, that those spin-induced oscillations 
are a direct conse quen ce of monopolc-quadrupole inter- 
actions m m, n, 11, . 
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All fully relativistic simulations presented here were 
performed using generalized harmonic coordinates. It 
would be very interesting to perform a similar study 
within the moving-puncture BSSN approach to inves- 
tigate whether our conclusions are applicable in other 
gauges. 
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